setwd("~/Documents/research/inProgress/philosophical/risk/")

## load and preprocess the data for the England Health Survey
source("replication_2017_09/rsr_ehsi_1loadprocess.R")

## load packages
library(ggplot2)
library(margins)
library(stargazer)

### main analyses: 

### Figures from the main text: 

## Figure 1: SR and subjective risk of job loss:

lj.sr2.0 = glm(sr.binary2 ~ losejob + eqinc.k + selfemp + female + age + I(age^2) , data = ehsi[which(ehsi$lookwk == 0 & ehsi$retired == 0 & ehsi$out_dis == 0 & ehsi$out_home == 0),],  family = binomial(link = "logit"))

#pd.lj2.0 = cplot(lj.sr2.0, "losejob")

par(mar = c(5,4,2,2)+0.2)
cplot(lj.sr2.0, "losejob", ylim = c(0, 0.6), xlim = c(-0.02, 1.02), xlab = "Subjective Probability of Job Loss", ylab = "Predicted Probability of Reporting Self-Respect", rug = F)

### Print command used for the figures in the paper (not run):
# dev.print(device = jpeg, width = 500, file = "subjR_sr_binary0.jpg")


## Figure 2: Objective risk: occupational unemployment rate

ue.sr2.0 = glm(sr.binary2 ~ occ_gen_uer + eqinc.k + selfemp + female + age + I(age^2) , data = ehsi[which(ehsi$lookwk == 0 & ehsi$retired == 0 & ehsi$out_dis == 0 & ehsi$out_home == 0),],  family = binomial(link = "logit"))

#ue.lj2.0 = cplot(ue.sr2.0, "occ_gen_uer")

par(mar = c(5,4,2,2)+0.2)
cplot(ue.sr2.0, "occ_gen_uer", ylim = c(0, 0.6), xlim = c(1.8,15.2), ylab = "Predicted Probability of Reporting Self-Respect", xlab = "Occupation-Gender Unemployment Rate", rug = F, xaxt = "n")
axis(1, labels = c("0.02", "0.04", "0.06", "0.08", "0.1", "0.12", "0.14"), at = c(2,4,6,8,10,12,14))
### Print command used for the figures in the paper (not run):
# dev.print(device = jpeg, width = 500, file = "uerR_sr_binary0.jpg")


##Figure 3 objective risk: skill specificity
sk.sr2.0 = glm(sr.binary2 ~ sksp.rel + eqinc.k + selfemp + female  + age + I(age^2) , data = ehsi[which(ehsi$lookwk == 0 & ehsi$retired == 0 & ehsi$out_dis == 0 & ehsi$out_home == 0),],  family = binomial(link = "logit"))

sk.lj2.0 = cplot(sk.sr2.0, "sksp.rel")

par(mar = c(5,4,2,2)+0.2)
cplot(sk.sr2.0, "sksp.rel", ylim = c(0, 0.6), xlim = c(0, 26), xlab = "Skill Specificity", ylab = "Predicted Probability of Reporting Self-Respect", rug = F)

### Print command used for the figures in the paper (not run):
# dev.print(device = jpeg, width = 500, file = "skspR_sr_binary0.jpg")

##### Figure 4: Expected probabilities for those who like risk, and those who don't like risk
lj.sr2b.risk = glm(sr.binary2 ~ losejob*riskseek + eqinc.k + selfemp + female  + marital + topqual2 + age + I(age^2) , data = ehsi[which(ehsi$lookwk == 0 & ehsi$retired == 0 & ehsi$out_dis == 0 & ehsi$out_home == 0),],  family = binomial(link = "logit"))

## 25th adn 75th percentile of risk orientation:
quantile(ehsi$riskseek, probs = c(0.25, 0.75), na.rm = T)

ld <- seq(0, 1, 0.025)

pr4 = predict(lj.sr2b.risk, 
	data.frame(losejob = ld,
   	riskseek = rep(4, length(ld)),
   	selfemp = rep(0, length(ld)), 
   	female = rep(0, length(ld)), 
   	age = rep(45, length(ld)),
   	eqinc.k = rep(mean(ehsi$eqinc.k, na.rm = T), length(ld))
   	), type = "response", se.fit = T)

pr1 = predict(lj.sr2b.risk, 
	data.frame(losejob = ld,
   	riskseek = rep(1, length(ld)),
   	selfemp = rep(0, length(ld)), 
   	female = rep(0, length(ld)), 
   	age = rep(45, length(ld)),
   	eqinc.k = rep(mean(ehsi$eqinc.k, na.rm = T), length(ld))
   	), type = "response", se.fit = T)

lower4 = pr4$fit - 1.96*pr4$se.fit
upper4 = pr4$fit + 1.96*pr4$se.fit

lower1 = pr1$fit - 1.96*pr1$se.fit
upper1 = pr1$fit + 1.96*pr1$se.fit


plot(c(0,1), c(0,0.7), type = "n", xlab = "Subjective Probability of Job Loss", ylab = "Predicted Probability of Reporting Self-Respect")

# se color in cplot: grDevices::gray(0.5, 0.5)

polygon(c(ld, rev(ld)), c(lower4, rev(upper4)), col = grDevices::gray(0.5, 0.5), border = grDevices::gray(0.5, 0.5))

polygon(c(ld, rev(ld)), c(lower1, rev(upper1)), col = grDevices::gray(0.5, 0.5), border = grDevices::gray(0.5, 0.5))


lines(ld, pr4$fit)
lines(ld, pr1$fit, lty = 2)

legend(0, 0.2, bty = "n", lty = c(1, 2),  legend = c("Like risk (75th percentile)", "Dislike risk (25th percentile)"))
legend(0.03, 0.09, bty = "n", fill = grDevices::gray(0.5, 0.5), legend = "95% confidence intervals", border = grDevices::gray(0.5, 0.5))

### Print command used for the figures in the paper (not run):
# dev.print(device = jpeg, width = 500, file = "subjR_likerisk_sr_DemandingBinary.jpg")
